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ALGORITHM  DEVELOPMENT  FOR  THE  TWO-FLUID  PLASMA 

MODEL 


AFOSR  Grant  No.  FA9550-05-1-0159 
U.  Shumlak 

Department  of  Aeronautics  and  Astronautics 
Aerospace  &:  Energetics  Research  Program 
University  of  Washington 


Abstract 

A  new  algorithm  is  developed  based  on  the  two-fluid  plasma  model  that 
is  more  physically  accurate  and  capable  than  MHD  ( magnet, ohydrodynamic) 
models.  The  algorithm  uses  high-order  spatial  and  temporal  accuracy  to  simu¬ 
late  time-dependent,  three-dimensional  plasma  phenomena.  High-order  spatial 
accuracy  is  accomplished  using  a  discontinuous  Galcrkin  finite  clement  method 
that  has  provided  up  to  lGt!l  order  accuracy.  The  temporal  evolution  is  ad¬ 
vanced  using  a  3r(i  order  Runge-Kutta  method.  The  numerical  fluxes  are  cal¬ 
culated  using  an  approximate  Riemann  solver  based  oil  the  two-fluid  plasma 
model.  The  source  terms  of  the  two-fluid  plasma  model  couple  the  electron  and 
ion  fluids  to  the  electromagnetic  fields.  The  simultaneous  solution  and  evolu¬ 
tion  must  be  tightly  coupled  to  prevent  unstable  numerical  oscillations.  Elec¬ 
tromagnetic  fields  arc  solved  by  both  formulating  Maxwell  s  equations  as  per¬ 
fectly  hyperbolic  equations  and  by  using  a  mixed  potential  formulation  which 
automatically  satisfies  the  divergence  constraint  relations.  Asymptotic  approx¬ 
imations  are  individually  applied  to  the  two-fluid  plasma  model  to  approach  the 
Hall-MHD  plasma  model.  An  improved  method  of  plasma  simulation  is  found 
by  using  the  two-fluid  plasma  model  with  an  artificially  increased  electron  to 
ion  mass  ratio  and  decreased  speed  of  light.  Multiscale  effects  arc  discovered 
in  cur rent-carrying  plasma  where  small-scale  electron  instabilities  lead  to  ion 
shocks  that  produce  large-scale  disruptions  on  the  plasma. 


1  Executive  Summary 

This  project  represents  a  three  year  effort  to  develop  a  new  algorithm  for  plasma 
simulations  based  on  the  two- fluid  plasma  model.  Current  plasma  simulation 
algorithms  capable  of  complex  geometries  are  based  on  the  MHD  (rnagneto- 
hydrodynamic)  model.  The  derivation  of  the  MHD  model  involves  several 
assumptions  that  severely  limit  its  applicability.  The  two-fluid  model  only 
assumes  local  thermodynamic  equilibrium  and  is,  therefore,  more  physically 
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accurate  and  capable  than  MHD  models.  The  two-fluid  model  is  formulated 
in  conservation  form.  An  approximate  Riemann  solver  is  developed  for  the 
two-fluid  plasma  model  to  compute  the  fluxes  in  a  stable  and  accurate  man¬ 
ner.  Several  methods  are  investigated  to  solve  the  electromagnetic  field  model, 
which  includes  the  source  terms  and  divergence  constraints.  These  methods 
include  a  characteristic-based  algorithm,  a  perfectly  hyperbolic  modification, 
and  a  mixed  potential  formulation.  The  two  plasma  fluids  and  the  electro¬ 
magnetic  fields  communicate  through  the  source  terms.  The  fluid  momentum 
and  energy  equations  have  source  terms  that  depend  on  E  and  B  The  elec¬ 
tromagnetic  equations  have  source  terms  that  depend  on  v*  and  ve  (Ampere’s 
law)  and  and  ne  (Gauss’s  equation).  Accurately  coupling  the  source  terms 
is  important  both  for  numerical  stability  and  for  modeling  plasmas  where  large 
equilibrium  forces  exist. 

An  algorithm  is  developed  for  the  complete  two-fluid  plasma  model  ini¬ 
tially  in  one  dimension.  The  algorithm  uses  a  Roc-type  approximate  Riemann 
solver  [1  to  discretize  the  hyperbolic  fluxes  of  the  fluid  model  and  an  upwind 
characteristic-based  solver  for  the  electromagnetic  fields.  Simulations  from  the 
resulting  finite  volume  algorithm  arc  benchmarked  against  known  analytical 
results.  Furthermore,  the  algorithm  is  applied  to  the  electromagnetic  plasma 
shock  problem  to  reveal  the  transition  from  gas  dynamic  shock  to  MHD  shock. 
The  results  arc  analyzed  to  reveal  the  fast  plasma  waves  that  arc  captured  in 
the  two-fluid  plasma  model.  [2 

A  high-order  algorithm  is  developed  that  uses  a  discontinuous  Galerkin,  fi¬ 
nite  clement,  method  for  the  spatial  representation  and  a  TVD  Rungc-Kutta 
method  for  the  time  advance.  [3]  Solutions  arc  found  with  up  to  16th  order  spa¬ 
tial  accuracy  and  3ni  order  temporal  accuracy.  The  two-fluid  plasma  algorithm 
is  used  to  model  niultiscale  physics  of  current-carrying  plasmas,  such  as  the 
Z-pinch  [4]  and  the  field  reversed  configuration  (FRC)  [5].  These  plasma  con¬ 
figurations  balance  large  equilibrium  forces  between  the  plasma  pressure  and 
the  electromagnetic  pressure.  The  high-order  algorithm  is  seen  to  significantly 
improve  the  ability  to  maintain  equilibrium  with  no  artificial  decay. 

The  divergence  constraints  of  Maxwell’s  equations  can  be  difficult  to  satisfy 
with  the  presence  of  current  and  charge  sources  on  an  arbitrary  computational 
grid.  The  divergence  constraints  are  satisfied  by  reformulating  Maxwell’s  equa¬ 
tions  to  include  correction  potentials.  The  approach  involves  coupling  the  di¬ 
vergence  constraint  equations  with  the  time-dependent  field  equations  to  form 
a  perfectly  hyperbolic  equation  set.  [6]  An  alternative  formulation  of  Maxwell’s 
equations  using  mixed  potential  is  also  implemented.  The  mixed  potential  for¬ 
mulation  automatically  satisfies  the  divergence  constraints. 

This  project  was  performed  by  Prof.  Uri  Shumlak  and  graduate  students 
Ammar  Hakim,  Robert  Lilly,  John  Loverich,  Bhuvana  Srinivasan,  and  Andrec 
Susanto.  This  project  resulted  in  doctoral  dissertations  and  master  thesis: 

•  John  Loverich,  “A  Discontinuous  Galerkin  Method  for  the  Two-Fluid 
Plasma  System  and  Its  Application  to  the  Z-Pinch”,  Ph.D.  2005. 
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•  Ammar  Hakim,  “High  Resolution  Wave  Propagation  Schemes  for  Two- 
Fluid  Plasma  Simulations”,  Ph.D.  2006. 

•  Bhuvana  Srinivasan,  “A  Comparison  between  the  Discontinuous  Galerkin 
Algorithm  and  the  High  Resolution  Wave  Propagation  Algorithm  for  the 
Full  Two-Fluid  Plasma  Model”,  M.S.  2005. 

These  dissertations  and  theses  can  be  obtained  from  the  University  of  Washing¬ 
ton  library  system  or  from  the  project  website,  http:// www.aa.washington.edu/cfdlab/. 
Archival  journal  and  conference  papers  were  published  reporting  on  the  work 
from  this  project: 

•  J.  Lovcrich  and  U.  Sliumlak,  “A  discontinuous  Galerkm  method  for  the 
full  two-fluid  plasma  model",  Computer  Physics  Communications  169 
251-255  (2005). 

•  A.  Hakim,  J.  Lovcrich,  and  U.  Sliumlak,  High  resolution  wave  propaga¬ 
tion  scheme  for  ideal  two-fluid  plasma  equations”,  Journal  of  Computa¬ 
tional  Physics  219  (1),  418-442  (2006). 

•  ,J.  Lovcrich  and  U.  Sliumlak,  “Non-linear  two-fluid  study  of  m— 0  sausage 
instabilities  in  an  axisymmctric  Z-pinclf  ,  Physics  of  Plasmas  13,  082310 
(2006). 

•  A.  Hakim  and  U.  Sliumlak,  “Two-fluid  physics  and  field-reversed  config¬ 
urations”,  Physics  of  Plasmas  14,  055911  (2007). 


2  Project  Results 


Plasmas  are  essential  to  many  technologies  that  are  important  to  the  Air  Force, 
some  of  which  have  dual-use  potential.  These  applications  include  portable 
pulsed  power  systems,  high  power  microwave  devices,  drag  reduction  for  hyper¬ 
sonic  vehicles,  advanced  plasma  thrusters  for  space  propulsion,  nuclear  weapons 
effects  simulations,  radiation  production  for  counter  proliferation,  and  fusion 
for  power  generation.  In  general,  plasmas  fall  into  a  density  regime  where  they 
exhibit  both  collective  (fluid)  behavior  and  individual  (particle)  behavior.  The 
intermediate  regime  complicates  the  computational  modeling  of  plasmas. 


2.1  Plasma  Models  —  Kinetic,  PIC,  MHD,  Two- 
Fluid 


Plasmas  may  be  most  accurately  modeled  using  kinetic  theory.  The  plasma  is 
described  by  distribution  functions  in  physical  space,  velocity  space,  and  time, 
/(x,v,£).  The  evolution  of  the  plasma  is  then  modeled  by  the  Boltzmann 
equation. 


Ofa 


dt 


(1) 
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for  each  plasma  species  a  =  ions,  electrons.  The  Boltzmann  equation  cou¬ 
pled  with  Maxwell’s  equations  for  electromagnetic  fields  completely  describe 
the  plasma  dynamics.  [7  9]  However,  the  Boltzmann  equation  is  seven  dimen¬ 
sional.  As  a  consequence  of  the  large  dimensionality  plasma  simulations  using 
the  Boltzmann  equation  arc  only  used  in  very  limited  applications  with  narrow 
distributions,  small  spatial  extent,  and  short  time  durations.  [10,11]  The  seven 
dimensional  space  is  further  exacerbated  by  the  high  velocity  space  that  is  un¬ 
used  except  for  tail  of  the  distribution  or  energetic  beams.  Boundary  conditions 
are  difficult  to  implement  in  kinetic  simulations. 

Particle  in  cell  (PIC)  plasma  model  apply  the  Boltzmann  equation  to  repre¬ 
sentative  supcrparticles  which  are  far  fewer  (107)  than  the  number  of  particles 
in  the  actual  plasma  (1020).  [12]  PIC  simulations  have  similar  limitations  as 
simulations  using  kinetic  theory  due  to  statistical  errors  caused  by  the  fewer 
supcrparticles.  Boundary  conditions  arc  also  difficult  to  implement  in  PIC 
simulations. 

The  other  end  of  the  spectrum  in  plasma  model  involves  taking  moments 
of  the  Boltzmann  equation  and  averaging  over  velocity  space  for  each  species 
which  implicitly  assumes  local  thermodynamic  equilibrium.  The  resulting  equa¬ 
tions  comprise  the  two-fluid  plasma  model.  The  two-fluid  equations  arc  then 
combined  to  form  the  MHD  model  [13]  However,  in  the  process  several  ap¬ 
proximations  are  made  which  limit  the  applicability  of  the  MHD  model  to  low 
frequency  and  ignores  the  electron  mass  and  finite  Larmor  radius  effects . 

The  MHI)  model  treats  the  plasma  like  a  conducting  fluid  and  assigning 
macroscopic  parameters  to  describe  its  particle-like  interactions.  Plasma  simu¬ 
lation  algorithms  based  on  the  MHD  model  have  been  very  successful  in  model¬ 
ing  plasma  dynamics  and  other  phenomena.  Codes  such  as  MACH2  are  based 
on  arbitrary  Lagrangian/Eulcrian  formulations.  [14]  ALE  codes  arc  well  suited 
for  simulating  plasma  phenomena  involving  moving  interfaces.  [15]  However, 
ALE  codes  cannot  be  formulated  as  conservation  laws  and  lack  many  of  the 
inherent  conservative  properties.  The  MHD  model  has  been  successfully  im¬ 
plemented  in  conservative  form  to  simulate  realistic  three-dimensional  geome¬ 
tries.  [1G,  17] 

A  severe  limitation  of  the  MHD  model  is  the  treatment  the  Hall  effect 
and  diamagnetic  terms.  These  terms  represent  the  separate  motions  of  the 
ions  and  electrons.  The  Hall  effect  and  diamagnetic  terms  also  account  for 
ion  current  and  the  finite  ion  Larmor  radius.  These  effects  are  important  in 
many  applications  such  as  electric  space  propulsion  thrusters:  Hall  thrusters, 
magnctoplasmadynamic  (MPD)  thrusters,  Lorcntz  force  thrusters.  The  Hall 
term  is  also  believed  to  be  important  to  electrode  effects  such  as  anode  and 
cathode  fall  which  greatly  affect  many  directly  coupled  plasma  applications. 
Furthermore,  the  Hall  and  diamagnetic  effects  may  be  important  for  hypersonic 
flow  applications.  [18] 

The  Hall  terms  can  be  difficult  to  stabilize  because  they  lead  to  the  whistler 
wave  branch  of  the  dispersion  relation.  The  phase  and  group  velocities  of  the 
whistler  wave  increase  with  frequency.  The  velocities  become  large  even  for 
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Figure  1:  Dispersion  relations  for  the  two-fluid  plasma  model  and  for  the  Hall-MHD 
plasma  model  that  results  when  asymptotic  approximations  are  applied  to  the  two- 
fluid  plasma  model  For  small  wave  numbers  and  low  frequencies  (right  plot),  the 
upper  branch  of  the  Hall-MHD  wave  follows  the  R  wave  of  the  two-fluid  model. 
However,  the  waves  diverge  and  the  Hall-MHD  wave  fails  to  follow  the  resonance 
at  the  cyclotron  frequency.  The  wave  speed  grows  without  bound.  Artificial  hyper¬ 
resistivity  is  required  to  damp  this  branch  of  the  Hall-MHD  wave. 


modest  values  of  the  Hall  parameter.  See  Fig.  1  for  the  dispersion  diagram 

A  semi-implicit  technique  has  been  applied  to  treat  the  Hall  term  in  a  Hall- 
MHD  model.  [19,  20]  After  the  hyperbolic  terms  of  the  MHD  equations  are 
advanced,  the  Hall  terms  arc  treated  independently.  The  conserved  variables 
arc  then  corrected.  The  procedure  can  be  computationally  intensive.  The  oper¬ 
ator  stencil  uses  5  points  in  the  sweep  direction  and  3  points  in  each  orthogonal 
direction.  The  complete  operator  stencil  is  45  points.  The  semi-implicit  method 
works  adequately  for  small  Hall  parameters,  but  becomes  unstable  or  slow  to 
converge  for  the  large  Hall  parameters  often  seen  in  applications. 

As  mentioned  above,  the  two-fluid  plasma  model  is  more  complete  than 
either  the  MHD  or  Hall-MHD  model.  The  two-fluid  plasma  model  resolves 
plasma  oscillations  and  speed  of  light  propagation.  However,  many  applications 
are  adequately  modeled  by  lower  frequency  dynamics.  Asymptotic  approxima¬ 
tions  (me  — »  0,  r  — ►  oc)  have  been  applied  to  the  two-fluid  plasma  model  to 
eliminate  the  high  frequency  waves  that  limit  the  maximum  numerical  time 
step.  Neglecting  electron  inertia  removes  the  limitation  due  to  the  electron 
plasma  and  cyclotron  frequencies.  Infinite  light  speed  removes  the  limitation 
due  to  light  transit  times.  The  asymptotic  approximations  reduce  the  two-fluid 
plasma  model  to  the  Hall-MHD  model  However,  applying  these  approxima¬ 
tions  fundamentally  changes  the  dispersion  relation,  as  evident  in  Fig.  1,  and 


6 


introduces  unphysical  wave  behavior.  Specifically,  the  phase  and  group  veloc¬ 
ities  of  a  Hall-MHD  wave  increase  without  bound  with  wave  number.  The 
large  wave  speeds  increases  the  stiffness  of  the  equation  system  making  accu¬ 
rate  numerical  solutions  difficult.  Furthermore,  the  maximum  wave  number  is 
usually  set  by  cither  the  computational  mesh  spacing  (kmax  ex  Ax)  or  by  an 
artificial  resistivity.  Rigorous  convergence  studies  are  difficult  with  the  simpler 
plasma  models  since  decreasing  Ax  leads  to  larger  kmax  and  shorter  wavelength 
phenomena. 


2.2  Two-Fluid  Plasma  Algorithm 


The  complexity  of  the  two-fluid  model  is  greater  the  MHD  model  but  signif¬ 
icantly  less  than  the  kinetic  model.  In  this  project  a  new  algorithm  is  devel¬ 
oped  that  solves  the  two-fluid  plasma  model  using  an  approximate  Ricmann 
solver.  ^2]  The  method  tracks  the  wave  propagation  across  the  domain  based 
on  conservation  laws. 

The  two-fluid  plasma  model  captures  the  separate  motion  of  the  electrons 
and  ions  without  the  added  complexity  of  the  kinetic  model.  The  two-fluid 
model  is  derived  by  taking  moments  of  the  Boltzmann  equation  for  each  species. 
The  process  of  taking  moments  eliminates  the  velocity  space  and  yields  rep¬ 
resentative  fluid  variables  (density,  momentum,  energy)  for  each  species.  The 
only  approximation  made  is  local  thermodynamic  equilibrium  within  each  fluid 
and  is,  therefore,  a  generalization  of  the  MHD  model.  The  fundamental  vari¬ 
ables  are  generated  by  taking  moments  of  the  distribution  function. 

The  evolution  of  the  particle  density  of  the  ions  and  electrons  is  expressed  by 
continuity  equations.  The  equations  are  the  zeroth  moment  of  the  Boltzmann 
equation. 


(2) 


(3) 


where  /i*,  ne  are  the  ion  and  electron  number  densities  and  the  particle  fluxes 
are  defined  by  the  partial  current  densities  jz  =  qn{vl  and  je  =  cneye  in  terms 
of  the  charges  q ,  e  and  fluid  velocities  Vj,  ve  for  the  ion  and  electron  fluids. 

The  first  moment  of  the  Boltzmann  equation  yields  momentum  equations 
for  each  species.  The  momentum  equations  arc  written  in  divergence  form. 
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where  E  and  B  are  the  electric  and  magnetic  fields,  pi  and  pe  arc  the  ion  and 
electron  partial  pressures,  and  R is  the  electron  to  ion  momentum  transfer 
vector. 
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The  second  moment  of  the  Boltzmann  equation  yields  energy  equations  for 
each  species  which  arc  expressed  in  divergence  form  for  the  total  energy. 


t +v 


(*  +  Pi)  — 
eni 


=  ji-  E  + 


Ref 


eri; 


(6) 


dfe 

dt 
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where  the  total  energy  is  defined  by 


£i 


-y; Pi  +  2nimiVi 


(7) 
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and 


^  = 


- 7Pe  +  -ne77leVg 

7  —  1  z 


(9) 


where  7  is  the  ratio  of  specific  heats  and  Ti,  Te  arc  the  ion  and  electron  tem¬ 
peratures.  A11  adiabatic  equation  of  state  is  assumed.  The  temperatures  are 
related  to  the  partial  pressures  by  pa  —  naTa  for  a  —  {z,  e}. 

The  equations  that  govern  the  ion  and  electron  fluids  arc  rewritten  in  com¬ 
pact,  conservative  form. 

^fVF  =  S  (10) 

where  Q  is  the  vector  of  conserved  fluid  variables,  F  is  the  tensor  of  hyperbolic 
fluxes  (Fx  +  Gy  +  Hz),  and  S  is  the  vector  containing  the  source  terms.  The 
vector  of  conserved  variables  is 


Q  =  [» 


'i  •,  1  Ji, t  1  Jiy  :  Jiz  1  Jex  i  Jay  i  Jez  i  i 
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(in 


for  the  number  densities,  electrical  current  densities,  and  energy  densities.  The 
flux  Jacobian  3F/3Q  for  the  two  fluid  equations  is  constructed  in  the  usual 
way.  The  characteristic  velocities  are  calculated  to  construct  the  approximate 
Riemann  fluxes. 

The  eigenvalues  of  the  flux  Jacobian  give  the  characteristic  velocities. 


The  eigenvalues  for  the  two  fluid  plasma  model  represent  the  combination  of 
the  drift  speeds  and  thermal  speeds  for  the  electrons  and  ions. 

The  electromagnetic  fields  influence  the  motion  of  the  plasma  fluid  through 
the  Lorcntz  force  which  is  contained  in  Eqs.  (4)  and  (5).  The  motion  of  the 
plasma  influences  the  evolution  of  the  electromagnetic,  fields  through  the  re¬ 
distribution  of  charge  density  and  current  density.  Maxwell’s  equation  govern 
the  evolution  of  the  electromagnetic  fields.  The  charge  density  qrii  —  ene  and 
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current  density  (j  —  jt  +  je)  are  calculated  directly  from  the  two-fluid  equations 
which  collide  the  electromagnetic  fields. 


<9B  „  „ 

—  =  -V  x  E 
dt 

(13) 

dE  .  . 

—  V  X  B//i0  -  (ji  +  Je) 

(14) 

eoV  •  E  =  qni  —  enc 

(15) 

V  B  =  0 

(10) 

In  multiple  dimensions  the  divergence  constraints  can  be  difficult  to  satisfy 
with  the  presence  of  current  and  charge  sources  on  an  arbitrary  computational 
grid.  The  divergence  constraints,  Eqs.  (15)  and  (16),  arc  satisfied  by  reformu¬ 
lating  Maxwell’s  equations  to  include  correction  potentials. 

The  approach  involves  coupling  the  divergence  constraint  equations  with  the 
time-dependent  field  equations  to  form  a  perfectly  hyperbolic  equation  set.  [6] 
The  field  equations  arc  expressed  as 


OB 

~dt 


4-  V  x  E  4-  7V  </>  =  0, 


(17) 


<9E  o  o  i 

--  -  c2V  x  B  +  xc2V<i>  = 

Ot  €0 

1  |  ^  E  _  (pij  -  ent. 

A  dt  <0 


_l_cty 

7c2  dt 


+  V  B  =  0, 


(18) 

(19) 

(20) 


where  0  and  are  the  electric  and  magnetic  correction  potentials  or.  more 
formally,  Lagrange  multipliers  which  vanish  at  the  domain  boundaries.  The 
method  more  accurately  predicts  the  propagation  speed  of  the  waves;  how¬ 
ever  the  method  can  overestimate  the  Lorcntz  force  caused  by  charge  sepa¬ 
ration.  The  implementation  illustrates  the  importance  of  tightly  coupling  the 
field  solver  to  the  fluid  solver. 

The  two-fluid  plasma  model  (including  the  electromagnetic  equations)  can 
also  be  expressed  in  divergence  form. 


8Q 

dt 


+  V  •  F  =  5 


(21) 


where  Q  is  t  he  vector  of  conserved  fluid  variables,  F  is  the  tensor  of  hyperbolic 
fluxes  ( F.r  4  Gy  4  Hz ),  and  S  is  the  vector  containing  the  source  terms. 

The  hyperbolic  fluxes  are  discretized  using  a  Roe-type  approximate  Rie- 
maiin  solver.  [1]  In  this  method  the  overall  solution  is  built  upon  the  solutions 
to  the  Riemann  problem  defined  by  the  discontinuous  jump  in  the  solution 
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at  cadi  cell  interface.  The  numerical  flux  at  the  cell  interfaces  is  written  in 
symmetric  form  as 

^+1/2  —  ^  (^i'+i  +  Fi)  -  2  ^  “  Qi)  1^*1  r*  C-2) 

where  /•*  is  the  ktfl  right  eigenvector,  A ^  is  the  eigenvalue,  and  4  is  the 
ktfl  left  eigenvector,  evaluated  at  the  cell  interface  (i  T  1/2).  The  values  at 
the  cell  interface  arc  obtained  by  a  Roc  average  of  the  neighboring  cells.  The 
flux  calculated  as  above  is  normal  to  the  cell  interface  which  is  the  desired 
orientation  for  applying  the  divergence  theorem  in  a  finite  volume  method. 

The  eigenvalues  and  eigenvectors  of  the  system  flux  Jaeobians  arc  calcu¬ 
lated  and  properly  normalized  to  prevent  catastrophic  cancellation.  [21,22]  A 
one-dimensional  approximate  Ricmann  solver  is  developed  based  on  the  de¬ 
rived  conserved  fluxes.  [2]  Electromagnetic  forces  arc  exerted  on  the  plasma 
fluids  through  the  source  terms  and  the  fluid  motion  affects  the  fields  through 
the  source  terms,  as  shown  in  Eq.  (21).  The  hyperbolic  fluxes  are  computed 
accurately  by  the  approximate  Ricmann  solver. 

Coupling  the  source  terms  to  the  hyperbolic  fluxes  is  critical  for  accurate 
simulations.  We  have  developed  a  wave  propagation  algorithm  that  using  a 
Strang  splitting  method  for  the  source  terms.  [23]  The  hyperbolic  fluxes  are 
computed  with  the  approximate  Ricmann  solver.  Limiters  used  on  the  hyper¬ 
bolic  fluxes  and  Strang  splitting  result  in  a  second-order  accurate  algorithm. 
However,  in  equilibrium  situations  where  forces  from  electromagnetic  fields  bal¬ 
ance  fluid  pressure  or  convective  forces,  the  contributions  from  the  source  terms 
must  be  accurately  calculated  to  balance  the  divergence  of  the  hyperbolic  fluxes. 
Even  small  errors  lead  to  a  diffusive  behavior.  The  source  terms  of  Eq.  (21)  are 
large,  in  general,  which  makes  the  equation  stiff. 

2.3  High-Order  Algorithm  for  Multiscale  Physics 

An  nnsplit,  finite-element  algorithm  is  developed  that  can  model  the  entire  two- 
fluid  plasma  model,  including  the  source  terms,  with  a  high  spatial  accuracy.  [3] 
The  algorithm  uses  a  discontinuous  Galcrkin  spatial  representation  with  a  TVD 
Rungc-Kutta  time  advance.  [24  26]  The  discontinuous  Galcrkin  method  is  a 
finite  element  approach  that  allows  for  arbitrarily  high  order  basis  functions  to 
model  the  variation  of  the  system  variables.  Source  terms  arc  automatically 
coupled.  The  method  currently  uses  up  to  sixteenth-order  accurate  spatial 
representation  with  a  third-order  accurate  TVD  Rungc-Kutta  time  advance 
method.  The  algorithm  has  been  implemented  for  multiple  dimensions  and  on 
parallel  computer  architectures. 

The  conserved  variables  of  the  two-fluid  plasma  model  arc  modeled  with  a 
set.  of  basis  functions.  {p/t}.  The  governing  equations,  expressed  as  Eq.  (21). 
are  multiplied  by  each  basis  function  and  integrated  over  the  mesh  cell  volume 
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(23) 


H.  Ail  integral  equation  is  generated  for  each  basis  function. 

[  vu^-dV  +  <f  vhF  dS-  [  F  •  VvhdV  =  [  vkSdV 

Jn  Jon  Jn  Ji  2 

whore  the  divergence  theorem  has  been  applied  to  the  second  term.  The  volume 
and  surface  integrals  arc  replaced  with  Gaussian  quadrature.  The  flux  F  is 
computed  with  the  approximate  Riernann  solver  with  a  limiter  applied  directly 
to  the  conserved  variables  to  get  high  resolution.  Less  accurate  Lax  fluxes  also 
typically  produce  adequate  results  with  reduced  computation  The  source  terms 
arc  described  by  the  basis  functions  and  are,  therefore,  the  same  order  accurate 
as  the  solution  variables.  The  high-order  representation  of  the  solution  variables 
satisfies  the  accuracy  requirement  to  preserve  the  equilibrium  balance  between 
the  divergence  of  the  hyperbolic  fluxes  and  the  source  terms.  Furthermore,  the 
source  terms  arc  directly  included  in  the  time  advance  of  the  solution  variables, 
and  no  source  splitting  is  necessary. 

The  discontinuous  Galcrkin  algorithm  has  been  applied  to  the  electromag¬ 
netic  plasma  shock  demonstrating  the  transition  from  gas  dynamic  shocks  to 
the  MHD  shock  [21,22]  as  the  Larmor  radius  is  reduced.  Analysis  of  the  data 
shows  the  differences  caused  by  the  additional!  plasma  waves  that  are  captured 
in  the  two-fluid  model  and,  consequently,  in  the  algorithm  developed  here.  [2] 
It  also  illustrates  the  dispersive  nature  of  the  waves  which  makes  capturing  the 
effect  difficult  in  MHD  algorithms.  The  electromagnetic  plasma  shock  serves 
to  validate  the  algorithm  to  published  data  (MHD  limit)  and  analytical  results 
(gas  dynamic  limit).  The  algorithm  has  also  been  applied  to  study  collisionlcss 
reconnection  and  the  results  arc  compared  to  published  results  of  the  GEM  chal¬ 
lenge  (Geospacc  Environmental  Modeling  Magnetic  Reconnection  Challenge) 
problem.  [27]  The  problem  is  difficult  to  model  and  provides  a  rigorous  test  for 
the  algorithm  and  benchmarks  to  other  algorithms.  The  evolution  of  the  re¬ 
connected  magnetic  flux  compares  remarkably  well  with  the  published  data  [3] 
Additional  applications  arc  discussed  in  more  detail  below. 

The  electromagnetic  field  model  includes  divergence  constraint  relations, 
which  if  not  accurately  satisfied,  can  lead  to  nonphysical  solutions.  Special 
treatment  is  required  because  the  divergence  relations  ovcr-constrain  the  solu¬ 
tion.  Satisfying  the  divergence  constraint  relations  can  be  accomplished  using 
a  Hodge  projection,  which  requires  solving  elliptic  equations  over  the  entire 
spatial  domain,  or  by  adding  correction  potentials  to  form  perfectly  hyperbolic 
equations,  which  requires  solving  additional  hyperbolic  equations  to  sweep  the 
divergence  error  out  of  the  domain,  as  described  above.  Alternative  to  these 
approaches.  Maxwells  equations  can  be  formulated  using  scalar  and  vector  po¬ 
tentials  that  automatically  satisfy  the  divergence  constraint  relations.  The 
mixed  potential  formulation  is  expressed  as 

d2(j>  2  02A  2  . 

—  -  v20  =  qm  -  ene,  -  V2A  =  j,  +  je,  (24) 

assuming  a  Lorentz  gauge  condition.  The  mixed  potential  formulation  is  im¬ 
plemented  by  solving  the  hyperbolic  equations,  given  by  Eq.  (24),  as  a  set  of 
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Figure  2:  Evolution  of  an  axisymmetric,  two- fluid  Z- pinch  with  an  initial  small,  si¬ 
nusoidal  perturbation  with  four  periods  showing  the  ion  density  contours.  Only  the 
final  stages  of  the  evolution  are  shown;  t  =  25,  30.  35,  40,  and  45.  The  perturbation 
remains  small  until  late  in  time  when  the  mode  becomes  unstable.  The  bending  of 
the  mode  is  an  expected  two-fluid  result  caused  by  the  finite  electron  mass.  Spatial 
scales  are  expressed  in  ion  Larmor  radii.  Vu/a  =  7.5. 


first-order  hyperbolic  equations  by  defining  auxiliary  variables.  The  gauge  con¬ 
dition  then  becomes  an  algebraic  expression  in  terms  of  the  auxiliary  variables. 
A  related  approach  is  to  assume  electromagnetic  waves  propagate  instanta¬ 
neously,  c  — *  oo,  known  as  the  Darwin  approximation.  [28]  This  approach  has 
been  implemented  as  part  of  the  asymptotic  approximations  described  above. 

2.4  Applications 

2.4.1  Multiscale  Structures  in  a  Z-pinch 

The  algorithm  has  been  applied  to  study  hybrid  plasma  instabilities  in  Z-pinch 
geometries.  [4]  The  results  arc  applicable  to  Z-pinches  experimentally  studied 
at  UW  and  Sandia  National  Labs.  An  axisymmetric,  two-fluid  Z-pinch  equi¬ 
librium  is  initialized  with  periodic  boundaries  in  the  axial  direction.  A  1%, 
sinusoidal  perturbation  of  the  azimuthal  magnetic  field  is  applied  and  the  plas¬ 
mas  dynamical  response  is  followed  The  effect  of  two-fluid  physics  can  be  seen 
by  adjusting  the  normalized  ion  Larmor  radius,  r^/a,  where  a  is  the  pinch 
radius.  Figure  2  shows  the  evolution  of  the  ion  density  for  during  the  final 
stages  of  an  instability;  t  =  25,  30,  35,  40,  and  45.  A  four  period  perturba¬ 
tion  is  applied.  The  applied  mode  is  seen  to  grow  eventually  entering  into  the 
nonlinear  regime  and  finally  plasma  confinement  is  destroyed.  The  protruding 
plasma  bends  downward  due  to  an  expected  two-fluid  phenomena  —  the  finite 
electron  inertia  influences  the  ion  density. 

When  the  plasma  size  is  reduced,  rn/a  —  2.5,  the  plasma  evolution  changes 
dramatically  as  seen  in  Fig.  3.  The  ion  density  contours  are  shown  at  t  =  10 
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Figure  3:  Evolution  of  an  axisyinmetric,  two- fluid  Z-pinch  with  an  initial  small,  sinu¬ 
soidal  perturbation  with  four  periods  showing  the  ion  density  contours.  The  earlier 
stages  of  the  evolution  are  shown;  t  =  10  and  15.  The  plasma  has  rn/a  —  2.5.  The 
perturbation  at  t  =  10  looks  similar  to  Fig.  2.  However,  at  t  =  15  the  four-period 
perturbation  is  overtaken  by  a  shorter  wavelength  mode. 


and  15.  The  same  equilibrium  and  perturbation  have  been  applied,  and  the 
four-period  perturbation  is  evident  at  t  =  10.  As  the  plasma  evolves  a  shorter 
wavelength  mode  grows  and  quickly  overtakes  the  initialized  perturbation.  The 
mode  is  independent  of  the  initial  perturbation.  The  wavelength  is  set  by  the 
Larmor  radius.  The  simulation  in  Fig.  3  has  rn/L  =  10,  so  a  ten-period  mode 
develops.  Furthermore,  the  finite  electron  inertia  leads  to  the  formation  of 
shocks. 

Analyzing  the  results  reveals  the  instability  occurs  when  the  electron  drift 
speed  exceeds  the  ion  sound  speed.  The  instability  is  related  to  the  lower 
hybrid  drift  instability.  t29]  If  the  electron  drift  speed  is  kept  below  the  ion 
sound  speed,  results  similar  to  Fig.  2  occur  if  the  equilibrium  is  MHD  unstable. 
However,  if  the  electron  drift  speed  exceeds  the  ion  sound  speed,  a  mode  similar 
to  Fig.  3  occurs.  This  effect  is  not  seen  in  MHD  plasma  models. 

2.4.2  Lower  Hybrid  Drift  Instability  in  a  Planar  Plasma 

The  algorithm  has  been  applied  to  study  hybrid  plasma  instabilities  in  field  Re¬ 
versed  configuration  (FRC)  geometries.  [5]  FRCs  are  experimentally  studied  at 
AFRL  and  UW.  The  same  effect  is  observed  in  planar  current-carrying  plasma 
sheets,  which  are  neutrally  stable  in  the  MHD  model.  The  simpler  geome¬ 
try  and  stability  properties  better  isolate  the  physical  mechanism  and  allow  a 
more  thorough  investigation.  An  axisyinmetric,  two-fluid  equilibrium  is  initial¬ 
ized  with  periodic  boundaries  in  the  longitudinal  direction.  A  single-period, 
sinusoidal  perturbation  of  the  density  is  applied  and  the  plasmas  dynamical 
response  is  followed.  Results  are  shown  in  Fig.  4.  The  spatial  scales  arc  nor¬ 
malized  by  the  ion  Larmor  radius.  The  initial  perturbation  docs  not  grow. 
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Figure  4:  Electron  density  evolution  of  a  current  sheet  at  t  =  100,  200.  250.  Current 
is  in-plane  to  the  left  with  a  confining  magnetic  field  out-of-plane  above  and  below 
the  current  sheet.  The  initial  equilibrium  is  MHD  stable,  but  develops  a  lower  hybrid 
drift  instability  captured  by  two-fluid  effects.  By  t  =  200,  perturbations  with  a  scale 
length  of  the  ion  Larmor  radius  are  visible.  The  instability  is  fully  developed  by  t  — 
250. 


Instead,  a  perturbation  with  a  scale  length  of  the  ion  Larmor  radius  develops. 
The  ions  separate  slightly  from  the  electrons  creating  an  electric  field  that  im¬ 
pedes  the  current.  The  finite  electron  inertia  leads  to  the  formation  of  shocks 
that  give  the  mode  a  “fishbone”  character. 

As  in  the  case  of  the  Z-pinch,  the  instability  occurs  when  the  electron  drift 
speed  exceeds  the  ion  sound  speed,  and  if  the  electron  drift  speed  is  kept  below 
the  ion  sound  speed,  instabilities  do  not  develop.  When  the  electron  drift  speed 
exceeds  the  ion  sound  speed,  the  mode  shown  in  Fig.  4  develops  even  though 
the  equilibrium  is  MHD  stable. 

The  instability  simulated  and  identified  has  practical  implications.  The 
lower  hybrid  drift  instability  has  been  suspected  as  the  cause  of  the  “anomalous 
resistivity”  in  FRC  experiments,  particularly  those  using  rotating  magnetic 
field  current  drive.  [30]  Numerically  computed  effective  impedance,  shows  an 
“anomalous  resistivity”  that  agrees  with  the  experimental  observations  and 
leads  to  cross-field  transport. 

2.4.3  3D  Instabilities  in  a  Z-Pinch 

Gross  three-dimensional  instabilities  in  a  Z-pinch  have  been  studied.  The  Z- 
pinch  results  are  shown  in  Fig.  5.  The  Z-pinch  equilibrium  is  expected  to  be 
unstable  to  gross  MHD  modes,  such  as  the  sausage  mode  (center  plot)  and  kink 
mode  (right  plot).  However,  an  additional,  small-scale  instability  develops  on 
top  of  the  MHD  modes.  The  instability  is  related  to  the  lower  hybrid  drift 
instability.  The  small-scale  structure  shown  in  Fig.  5.  is  not  captured  in  the 
Hall-MHD  and  MHD  plasma  models. 
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Figure  5:  Electron  density  evolution  of  a  Z-piiich  showing  the  development  of  a  lower 
hybrid  drift  instability  superimposed  on  the  sausage  and  kink  MHD  modes. 

2.5  Conclusions 

A  motivation  for  implementing  the  high-order  discontinuous  Galcrkin  method  is 
to  accurately  capture  the  detailed  spatial  structure  of  plasma  dynamics  without 
necessitating  large  computational  grids.  We  have  investigated  this  ability  by 
comparing  solutions  from  the  discontinuous  Galcrkin  method  with  a  second- 
order  wave  propagation  method  applied  to  a  variety  of  hyperbolic  problems 
linear  advection,  electrostatic  ion  cyclotron  waves,  electromagnetic  waves,  and 
two-fluid  plasma  dynamics.  The  general  finding  indicates  that  for  applications 
with  a  single  characteristic  speed  or  speeds  with  a  limited  range,  the  low-order 
method  adequately  captures  the  solution  with  substantially  less  computational 
effort.  However,  for  applications  with  disparate  characteristic  speeds,  the  high- 
order  method  is  better  able  to  capture  the  solution  without  the  phase  errors 
that,  appear  in  the  low-order  method. 

The  two-fluid  plasma  model  resolves  plasma  oscillations  and  speed  of  light 
propagation.  However,  many  applications  arc  adequately  modeled  by  lower 
frequency  dynamics.  Asymptotic  approximations  ( me  — >  0,  c  — >  oo)  have  been 
applied  to  the  two-fluid  plasma  model  to  eliminate  the  high  frequency  waves 
that  limit  the  maximum  numerical  time  step.  Applying  these  approximations 
fundamentally  changes  the  dispersion  relation  and  introduces  unphysical  wave 
behavior.  Accurate  simulations  require  large  computational  efforts. 

More  accurate  and  less  computationally  intensive  simulations  arc  possible 
using  the  two-fluid  plasma  model  with  reduced  mass  ratio  and  light  speed.  The 
high  frequency  dynamics  captured  by  the  two-fluid  plasma  model  is  modified 
by  allowing  the  electron  mass  to  increase,  such  that  the  mass  ratio  rrii/m e  is 
smaller  than  the  physical  value  (1836  for  hydrogen  plasma).  We  have  produced 
accurate  simulations  with  mass  ratios  as  small  as  25.  Increasing  the  mass  ratio 
to  100  docs  not  significantly  alter  the  results.  Similarly,  the  speed  of  light 
only  needs  to  be  much  larger  than  the  next  fastest  characteristic,  typically  the 
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Alfven  speed.  The  ratio  c/vj\  is  approximately  1000  in  experimental  plasmas. 
However,  we  have  produced  accurate  simulations  with  values  as  small  as  10. 
These  reduced  values  provide  significant  computational  speed  increases  without 
a  significant  loss  of  accuracy. 
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